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Abstract 

^ Monte Carlo sampling of nonlinear state-space models is particularly difficult 

, c^, in circumstances where the transition density is not of closed form. This is com- 

mon in the physical and biological sciences, where real phenomena demand mod- 
^ elling that accurately reproduces nonlinearity, chaos, mass conservation, multiple 

potentials and other behaviours that do not necessarily yield convenient analyt- 
ic ical forms. In this work, we explore a new idea in the design of samplers for 
^ such models, collapsing out the state variables of a state-space model to leave only 
its noise terms. We then consider the design of proposal distributions over these 
^ noise terms, exploiting the independence and simple parametric forms prescribed 

to them in the prior structure of the model. 
^ The chosen vehicle for joint state and parameter estimation in these collapsed 

>■ state-space models is the particle marginal Metropolis-Hastings (PMMH) sampler, 

^ from the family of particle Markov chain Monte Carlo (PMCMC). An auxiliary 

particle filter is used in the inner loop, with proposal distributions over noise terms 
tuned using the unscented Kalman filter. We conduct a thorough empirical investi- 
gation of the PMMH sampler in this context, including a look at the improvement 
obtained in acceptance rates, convergence rates and variability of the likelihood 
estimator. Case studies from the domain of marine biogeochemistry are used to 
demonstrate the ideas. We believe that these cases, exhibiting mass conservation 
and mild chaotic behaviour, are particularly challenging applications on which to 
exercise the PMCMC methodology. 

Keywords: marine biogeochemistry, particle filter, particle Markov chain Monte Carlo, 
sequential Monte Carlo, unscented Kalman filter 
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1 Introduction 



The tuning of a proposal distribution to improve the performance properties of a 
Monte Carlo sampler requires both a prior density and likelihood function that are 
computationally tractable. In rejection, importance and Markov chain Monte Carlo 
(MCMC) sampling, computations appear of the form: 

likelihood x prior density 
proposal density 

In the case of a state-space model (detailed further in the prior density is: 

T 

pi^O:T,0) = p{^o\O)p{e)l[p{^t\^t~U0) , (2) 

t=l 



and is usually specified via explicit densities p(xo|0) and p{0), with the remainder im- 
plied by the transition density p(xt|xj_i, 0). If the transition density has a convenient 
closed form, such as for conditional Gaussian models, then the prior distribution is 
computationally tractable. This facilitates sampling using arbitrary proposal distribu- 
tions. Across the spectrum of potential applications for nonlinear state-space models, 
however, such cases may well be exceptional, not the common norm. If the transition 
distribution can be simulated, but has no closed-form probability density function [see 
e-g-EQUsHlll, then the prior distribution is computationally intractable. In such cases, 
one approach is to constrain the candidature of proposal distributions to those that 
cancel the prior density (or at least its intractable factors) in ([T]). Simulation from the 
prior must feature in the proposal to achieve this, as in the bootstrap particle filter [llj. 
An alternative is to unbiasedly estimate the transition density, as in the random-weight 
particle filter [8J. Both approaches can be quite limiting in a multivariate setting. 

Models without closed-form transition densities are common in the biological and 
physical sciences, where nonlinear functional forms of increasing complexity are used 
to capture various phenomena in ever-greater detail. Such models may feature bi- 
furcation, multi-well potentials or chaos, and stochasticity may enter by continuous, 
currency conserving or other elaborate means to render unwieldy transition densities. 

The first contribution of this work is to remedy the proposal problem by collapsing 
the state-space model to its noise terms: the independent, univariate terms of pre- 
scribed distribution that form the source of stochasticity (Figure [T]). This can be done 
for any model, or at least, pragmatically, must be doable if we intend to rely on se- 
quences of independent (pseudo-)random numbers. The effect is that the likelihood 
function subsumes the nonlinear complexity of the model, leaving a simple prior den- 
sity over noise terms. The transition density becomes trivial, noise terms being inde- 
pendent a priori, and one can begin to design proposal distributions over these. We 
should caution that the design of such proposals may be difficult, a strongly nonlin- 
ear mapping between noise terms and observations potentially producing complex 
posterior distributions over those noise terms, which may be difficult to approximate. 
Where a closed-form transition density is available, collapsing to noise terms is un- 
likely to be of benefit. Nevertheless, one might appreciate the ability to exert at least 
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Figure 1: Graphical models of (a) conventional state-space model, and (b) collapsed 
state-space model, where the conventional model is reduced to initial conditions, noise 
terms and observations. Parameters, ©, have been removed for clarity, but note that 
Ui:T -11 ©, while all other variables depend on them. 



some control, however limited, for a class of models where inference procedures are 
usually quite constrained. 

By way of context and completeness, an intractable likelihood function, rather than 
prior density, would lend itself to Approximate Bayesian Computation [|29||. The pro- 
posed collapse to noise terms differs from the separation of tractable and intractable 
components in the Rao-Blackwellisation of state-space models [5J, where it is typical 
to admit dependence of the tractable (linear) component on the intractable (nonlin- 
ear) component, but not vice-versa; in the collapsed state-space model described here, 
the intractable (nonlinear) component depends on the tractable (independent noise) 
component. 

In this work, the particle marginal Metropolis-Hastings (PMMH) sampler is em- 
ployed, a particular particle Markov chain Monte Carlo (PMCMC) algorithm 111. The 
unscented Kalman filter [UKF, |15l |39l is used to construct the proposal distributions 
via which improved performance of that PMMH sampler is sought. Two models from 
the domain of marine biogeochemistry are used as a study for this, previously treated 
with simpler PMMH approaches in Jones et al. [14J and Parslow et al. [25J. To tame 
the PMMH sampler, we must better understand its behaviour in the wild: quite some 
space is dedicated to this, and the insight and diagnostics developed form the second 
major contribution of the work. 

The collapsed state-space model is described in ^ PMMH, along with prelimi- 
naries on its behaviour, are described in ^ Particular proposal mechanisms within 
PMMH using the UKF are given in ^ Case studies feature in ^ with detailed dis- 
cussion and conclusions in ^ Appendices pick up the statistic used for assessing the 



likelihood estimator across space (^ A), efficiencies for the special case of models with 
only Gaussian noise (^, along with other non-essential implementation issues, and a 
description of the substantive model used in the second case study ([Q. 
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2 The collapsed state-space model 



For T time points, a sequence of observations yi, . . . , yr of random variables Yi, . . . , G 



y is assumed given, indicative of latent states Xi , . . . , G M and initial condi- 
The model is parameterised by static © G M^'^, with the hierarchical 



tion Xq G 
model: 



0) 



t=i 



JJp(Xt|Xi_i,0) 



t=l 



(3) 



We follow the convention that uppercase letters denote random variables, while low- 
ercase letters denote values of them. Bold typeface indicates vectors. For some ran- 
dom variable A, the notation A^.j denotes the sequence A^, A^+i, . . . , A^, and A]^ " the 
sequence of A^ samples, A|", A™"*"^, . . . , A". 

This standard state-space model formulation is depicted in Figure [l][a), along with 
the independent noise terms Ui, . . . , Uy G M^" that drive the transition. In the con- 
text of the models under study, it is expedient to collapse the latent state to these 
noise terms only, as in Figure [T](b), where the Xi ^- are eliminated as mere intermedi- 
ates in the mapping from {Xq, Ui^^, ©} to {Yi-^}. This permits the design of sam- 
plers over the random variables {Xq, Ui:t, ©}/ rather than the typical set {Xo;t, 0}- 
A closed-form transition density is assured by the independence prescribed to Ui:t in 
the model's prior structure. In this work, the Ui:t are always standard i.i.d. A/'(0, 1) 
variates, independent of parameters in the prior structure, but scaled and translated 
by them as required. Variates of any base distribution might be equally used, however. 
The new hierarchical model is 



P(yi:T,Ui:T,Xo,0) 



t=l 



t=i 



p{^o\o)p{e) 



(4) 



Simulation of the process model is described by the function / : {Ut, X(_i, ©} — > 
Xt. Recursive application of /(■) can be used to recover the state trajectory xi:^ from 
any sample {ui-^, xq, 6}. 

For simplicity, we assume that A^^., A^^ and A^^^ are constant across time. Note, how- 
ever, that the methods do not depend on this, or even on equispaced observations; 
they are equally applicable with changing state and observation sizes, such as under 
missing or sparse data. 



3 Inference methods for the collapsed state-space model 

The target (posterior) density is 7r(ui:T, xq, 0|yi:T)/ factorised as either: 

7r(ui;T,xo,0|yi:r) = 7ri(ui:r, Xo|0, yi:T)7r2(0|yi:T) • (5) 

or 

7r(ui:r,xo,0|yi:r) = vri(ui:T|xo, ^, yi:T)7r2(xo, 0|yi:r) ■ (6) 
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In either case, the first factor, 7ri(-), is targeted using the formulation of the auxiliary 



particle filter [APF, |26l, described in ^:3.1 The second factor, tt2{-), is targeted in an 



outer loop around the particle filter using Metropolis-Hastings [MH,'2T1|T2J. The par- 
ticle filter nested within MH defines particle marginal Metropolis-Hastings [PMMH, 
[11, described in §3.2[ 

The factorisation ^ is attractive over (jSj) in circumstances where prior informa- 
tion over initial conditions is scarce, and an importance sample over them would be 
degenerate. It is used for the derivations here; use of the alternative is straightforward. 

3.1 Auxiliary particle filtering 

For given {xq, 6}, estimation of 7ri(ui.t|xo, 6, yi,t) iriay proceed recursively as: 

p{ui,t\xo,0,yi;t) oc j9(yt|ui:t,xo,0)p(ui,t|xo,0,yi:t_i) (7) 
= p{yt\ui:t, xo, 0)p{ut\ui.,t^i, xo, 6, yi:t_i)p(ui:t_i|xo, 6, yi:t-i) (8) 
= p(yt|ui:i,xo,0)p(ut)p(ui;t_i|xo,0,yi:i_i) . (9) 

The key step, and motivation for operating over {Ui:t, Xq, ©}, is the third. An or- 
thodox derivation over {Xo:r, ©} yields the transition density p(xt|xt_i, 6), for which 
a closed form is not expected in the models under study. The collapsed state-space 
model instead produces the term p{ut), of closed form by definition. 

At time t, a sequential importance sampling of (j9]) proceeds with a draw u™^ ~ 
qt{ui,t)> weighted with 

m _ P(yt|Ubi> Xq, 6>)p(u^)p(u^,_i|xo, 6, Ji^t^i) 

Here, qt{iii,t) is some proposal distribution satisfying the usual assumptions for se- 
quential importance sampling [4J. 

The reduction to noise terms does not produce a method profoundly different to 
the conventional. The construction is simply that of the auxiliary particle filter [26J, 
reduced to {Ui:^, Xq, ©}. State samples x™ = f{u^,±^-^^,6) may be stored and re- 
cursively updated in place of complete trajectories of noise terms u™^, and no densi- 
ties that would require a Jacobian correction for change-of -variable are ever explicitly 
evaluated. The only change is to enable arbitrary proposal distributions in the space 
of noise terms. 

A fully adapted ||^l28ll proposal would take the form 

gt(Ul:t) = gi(ui:f_i|xo,0,yi:t)g2(Ut|Ui:t„i,Xo,0,yi:t) . (11) 

It is not expected that full adaptation can be obtained in general, but reduction to 
noise terms does facilitate partial adaptation - various approximate substitutes for 
gi(-) and g2(')- Sampling from gi(-) typically involves performing some auxiliary look- 
ahead procedure to time t IITSlI , yielding weights from which a look-behind resam- 
ple at time t — 1 is performed to produce samples. The resampled particles are then 
propagated forward, with proposals drawn from g2(-)- refer to the weights w^'*^ 
as stage-one weights, and the final weights wl'-^^ as stage-two weights. The concepts 
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of look-ahead and look-behind are purposely distinguished here to avoid confusion 
later: some methods will feature a look-ahead using a UKF to construct q'2(-)/ without 
then using a look-behind to construct qi{-). 

Pitt [27J considers estimation of the likelihood p(yi;r|^) with the APF, equally ap- 
plicable to p(yi:T|xo, 0)- The estimator uses both stage-one and stage-two weights, and 
is given by 

p(yi:Tixo,e)-n|]^E<|{E^r| • (12) 

t=l [_ m=l J Lm=l J 

Proofs of unbiasedness are given in Del Moral BH and Pitt et al. |I28II . 



3.2 The particle marginal Metropolis-Hastings sampler 

At a high-level, the PMMH sampler combines a MH sampling over {Xq, ©} with a 
sequential Monte Carlo (SMC) sampling over {Xi-^^lXo, ©}, or {Ui:r|Xo, ©} in the 
collapsed case. In the outer loop, a proposed move {xq, 6} — )■ {xq, 6'} ~ ^(xq, 0'|xo, 6) 
is accepted with probability 



mm 



p{y,.,TK,e')p{^'„6')q{^o,6\^'„0') 



(13) 



where the likelihoods are estimated by SMC targeting 7ri(-) in ([6]); in this work by an 
APF as described above. 

Rigorously, 7ri(-) is a marginal of a distribution over an extended space in which 
the PMCMC algorithm operates, a space which includes variables associated with the 
resampling mechanism of the particle filter [c.f . Equation 22 of [U: 



l:M 
T 5 



1:M 
'1:T 



xo,e)=n 



t=i 



M 
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l:M 



W- 
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u 



1 ) 



u. 



-1 

-1 5 



xo,6l) 



(14) 



Here, a]^ is the index of the particle at time t — 1 which is the ancestor of particle 
m at time t, as selected by the resampling algorithm. By recursively tracing these 
indices backward to time s, where s < t, the notation 6™ denotes the index of the 
particle at time s which is the ancestor of particle m at time t. For some specific 
m ~ Multinomial(w^ ), a sample of 7ri(-) is given by the marginal {u^^ , • • • , }, 
the precise form of which is not required. 

PMMH chains are disposed to stammering if the likelihood estimator ( [T2| > has high 
variability. A chain moving into a particular state on the basis of an unusually large 
likelihood estimate will tend to remain there for a prolonged period, resisting even 
very local proposals; it will finally concede to a new state after a perhaps inordinate 
amount of time. The source of variability is both sampling and resampling error in 
the particle filter ||27l[T6ll . The methods proposed in this work target a reduction of the 
former for any given number of particles. 

In Pitt et al. [28J, one measure of performance used to compare algorithms is the 
standard deviation of repeated log-likelihood estimates at the true parameter value. 
The truth is not, of course, known for real observational data, but nevertheless it 
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is helpful to perform a few pilot runs at carefully selected 6 and tune M until the 
statistic is tolerable. While a satisfying configuration might be found for particular 
regions, rarely is one achievable across the full expanse of the prior: the estimator can 
be strongly heteroskedastic, and not independent of the parameters themselves! The 
problem is particularly acute when the process model itself is used as part of the pro- 
posal within the APR Small values of those parameters which prescribe process noise 
variance, for example, may lead to narrow proposal distributions that amplify weight 
variance, and so the variability of likelihood estimates. The mixing properties of the 
model may also be affected by, for example, gradient- and decay-related parameters. 

To explore this behaviour, we approximate, for any sample {xq, 6}, the long-term 
acceptance rate of repeated proposals at that point. We refer to this as the conditional 
acceptance rate (CAR) at {xq, 0}, with details given in ^|Aj By computing this statistic 
at many points across the space of parameters and initial conditions, it is possible to 
explore the heteroskedasticity of the likelihood estimator. While related to the alter- 
native metric of standard deviation Il28l , the CAR is more directly interpretable as to 
the impact of variability on the acceptance rate of a chain, and also accommodates any 
asymmetry in that variability. We prefer it for these reasons. 



4 Proposal mechanisms for the collapsed state-space model 

The combination of the auxiliary particle filter and collapsed state-space model per- 
mits the design of arbitrary proposal distributions, even for those models which would 
otherwise have no closed-form transition density. Three proposal schemes are consid- 
ered in this work, each with and without a look-behind, to give six methods in total. 
The first scheme is the ordinary bootstrap particle filter ||TT|. The remaining two use 
the UKF in different ways, one similar to the existing unscented particle filter [37], al- 
beit in noise space. All six methods are detailed in this section; consult Table [l]for the 
quintessential. 

4.1 Bootstrap particle filter 



The simplest approach, that of the bootstrap filter ||TT|. sets 

= p(ui:t_i|xo,0,yi:f_l)p(Ut) . (15) 

This is straightforward to apply, and indeed does not require a collapse to noise terms 



for tractability. Note that both factors of ( [15l > cancel in the weight computation of ( 10 1, 
so that particles are merely weighted by likelihoods. We denote this method PFO. 

Lacking p(yf |ui.(_i, xq, 0), an analytical look-ahead is not forthcoming. A deter- 
ministic single-point pilot look-ahead, simulating with Uj = 0, can offer improvement 
in some circumstances [18J. We denote such a method PFl: 

gt(Ul:t) = p(ui:f_i|Ut = O,Xo,0,yi:t)p(Ui) (16) 
OC p(yi|Uj = O,Ui:t_i,Xo,0)p(Ui:t_i|xo,0,yi:t_l)p(Ui) . (17) 
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4.2 Marginal unscented particle filter 



We next attempt to draw on analytical approximations to inform the proposal distri- 
bution. The particular focus is on the UKF, which, for modest state sizes, tends to 
outperform other approximate nonlinear Kalman filtering approaches p9|, such as 
the extended and ensemble BHIZ!] variants. The UKF approximates the time-marginals 
p{ut, Xj |xo, 6, yi,t) using a Gaussian distribution. We denote the approximation p^(ut, |xo, 6, yi,t] 
henceforth. At each time, 2(Nu + A^^: + A^^) + 1 number of cr-points are crafted about 
the mean of the Gaussian distribution, propagated through the process model and 
specifically weighted to compute pAr(ut, xt, y^lxo, 0, yi:t-i). Conditioning this on the 
actual observed value yt delivers the approximate time marginal pxi'^t, Xf|xo, 6, yi:t). 
See Julier and Uhlmann [151 arid Wan and van der Merwe [39J for details. 

The first UKF-based approach adopted is to use the marginal UKF approximations 
PA^(ut|xo, 0, yi;t) at each time as a common proposal for each particle. We call this the 
marginal unscented particle filter (MUFF), and without look-ahead, specifically MUPFO: 



gt(ui:i) = p(ui;t_i|xo,0,yi:f_i)pAr(ut|xo,0,y 



l:t 



(18) 



We can again consider a deterministic single-point pilot look-ahead, although this time 
might assume that fixing to /ij, the mean estimate obtained from the UKF, will 
provide a better result than fixing to 0. We denote this MUPFl: 

gi(ui:t) = p(ui:(_i|ut = /ii,xo,0,yi:t)pA^(ut|xo,0,yi:t) (19) 
oc p(yt|uf = /ii,ui;t_i,xo,0)p(ui:t_i|xo,0,yi:t-i)pAr(ut|xo,0,yi:t) . (20) 



Note that (18l and (20l use the same time-marginal pAr(ut|xo, 0,yi:t) for each par- 
ticle, not the conditional for the mth particle, pJj-{\i"'-\vL^j^_i, xq, 0, yi;t), which may be 
preferred conceptually (and is the basis for the next scheme, see ^ 4.3[ >. This is a compu- 
tational concession, permitting the UKF to be run asynchronously to the particle filter 



(see ^B.l I, but is not without justification for fast-mixing models regardless. The over- 
head of the method is slight, with only 2{Nu + + Ny) + 1 additional propagations 
for each particle filter. This is likely small compared to M, the number of propagations 
of the particle filter, which typically scales exponentially with dimension. 

In the case where Ui;^ ~ i.i.d. A/'(0, 1), the also-Gaussian marginals of the UKF 



afford some efficiencies in the weight calculation after slotting (18 1 or (20 1 into (10 1. 
These are given in ^B.2 This also suggests that, even in the presence of a closed-form 
transition density over state variables, a collapse of the state-space model to noise 
terms may still be of benefit, computationally, to speed up weight calculations. 



4.3 Conditional unscented particle filter 

Separate UKF look-aheads for each particle may be preferred to a marginal lookahead, 
particularly for slower-mixing models. We refer to this as the conditional unscented 
particle filter (CUFF). It is not dissimilar to the unscented particle filter 1 371 , besides the 
collapse to noise terms. Without look-behind, the CUPFO method uses: 

qt{Vil:t) = p(ui:t_i|xo, 0, y i:i_i)pAr(Ut | Ui:t_i , Xq, 0, yi.t) . (21) 
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Note that the proposals are now different for each particle, as, after drawing from 
the first component, a UKF must be run to construct the second component for sam- 
pling. The proposal is a substantially more expensive operation as a result. Each UKF 
requires 2{Nu + Ny) + 1 number of a-points. While this is fewer than the MUFF (as 
the starting state for each is the single point x™ j^), the total number of propagations 
required is nonetheless substantially more, at 2M{Nu + Ny + 1). Furthermore, the 
UKFs must be synchronous to the PF, eliminating the opportunity for asynchronous 
execution on multiple devices as described for the MUFF methods in ^B.l As conso- 



lation, weights for the look-behind are essentially free, giving a CUPFl method that 
potentially improves substantially on CUPFO at little additional cost: 

gt(ui:i) = pAr(ui;t_i|xo,0,yi:t)pAr(ut|ui;t_i,xo,0,yi:t) (22) 

OC PA^(yt|ui:t_i, Xo, 0)p(ui:t_i|xo, 0, y )pAr (Uf | Ui:t_iXo, 0, Yl-t) ■ (23) 

For higher-dimensional models the complexity may be dominated by the computation 
of the Schur complement required in the Kalman update equations. This is simplified 
by the independence of noise terms, however, with the particular case for Gaussian 
noise treated in §B.2 



5 Case studies in marine biogeochemistry 

The proposed methods are assessed empirically on two models in the domain of ma- 
rine biogeochemistry. Comparisons include conditional acceptance rates (CAR) and 
the spatial dependence of these, overall acceptance rates and Markov chain conver- 
gence rates. 

Methods are assessed in two configurations (see Table |2|): particle-matched, where 
the number of particles, M, is the same for all methods, and compute-matched, where M 
is set for all methods so as to roughly equate the number of propagations performed, 
where these include look-ahead pilots and UKF a-points. The compute-matched case 
is meant to control, to some extent, for the particular software implementation of each 
algorithm. Rather than matching wallclock execution time, it instead recognises that 
propagations typically, and almost universally, dominate the profile of any particular 
implementation (80-90% of execution time in the authors' experience). 

Simulated data is used, generated from the case study models themselves. This 
has a number of advantages over real observational data: execution time can be man- 
aged to facilitate the many runs required for some diagnostics, the model is perfect 
in the generative sense, capturing all (and only) those processes influencing observa- 
tions, and a known ground truth for all latent and observed variables is available for 
validation of the methods. Note immediately, however, that given the simulated data 
set has noisy and finite number of observations, a hypothetically perfect Bayesian in- 
ference will not recover the ground truth exactly: not only is the ground truth not the 
maximum likelihood estimate (MLE), but the prior may bias the result away from this. 

An observational data set is available, specifically from Ocean Station Papa (OSP), 
to which the second model has previously been fit using similar, but simpler, meth- 
ods [25 J. While the methods presented here can be used on this, the story for this 
data set is as yet incomplete, particularly around rigorous handling of its sparsity 



10 









Put tic le-ma tched 




Cotnpute-tnatched 


v^ase study 


Metnou 


M 


Acceptance rate 


M 


Acceptance rate 


PZ 


PFO 


64 


.182 (.003) 


384 


.245 (.002) 




PFl 


64 


.189 (.003) 


192 


.233 (.002) 




MUl rU 


d4 


.zUo (.UUJj 


J/fa 


.ZDi {.(JUZ) 




iViUi r i 






io4 


.z4o {.(jVZ) 




(^Ul rU 


d4 


.zUy {.(JUo) 


d4 


.zUy [.UUo) 




(^Ul ri 


o4 


.zi4 [.(jUo) 


o4 


.zi4 (.UUJ; 


NPZD 


PFO 


64 


.142 (.003) 


1536 


.276 (.002) 




PFl 


64 


.139 (.006) 


768 


.254 (.003) 




MUPFO 


64 


.165 (.003) 


1504 


.278 (.002) 




MUPFl 


64 


.167 (.007) 


736 


.263 (.003) 




CUPFO 


64 


.169 (.003) 


64 


.169 (.003) 




CUPFl 


64 


.177 (.005) 


64 


.177 (.005) 



Table 2: Configuration of methods for case studies, with mean (and standard devia- 
tion) of resulting acceptance rates across 256 chains. 

(up to seven weeks between observations, compared to daily for the simulated sets 
here). Until that story is complete, we have chosen not to draw on this data for testing 
purposes, but pick up on the point in discussion (Q. Besides this issue of sparsity, 
however, we believe that the second case study is quite representative of real-world 
cases, especially in its nonlinear functional forms. 

5.1 PZ model 

The first model considered is a variant of the Lotka-Volterra differential system Ill9ll38ll , 
specifically over the predator-prey relationship of zooplankton and phytoplankton in 
a marine environment. Previously treated with a PMMH-style sampler [14], the intent 
here is to plumb deeper into the behaviour of the algorithm, the presence of just two 
parameters providing an ideal opportunity to visualise dependence of CAR on ©. This 
PZ (phytoplankton and zooplankton) model modifies the classic Lotka-Volterra with 
the addition of a quadratic mortality term for zooplankton and stochastic growth term 
for phytoplankton. The stochasticity admits varying growth rates in phytoplankton 
without explicitly modelling contributory factors such as light and temperature, and 
thus exemplifies how such uncertainties can be treated by the introduction of stochas- 
ticity into an otherwise deterministic model [|T4ll25l . 

The state of the model is given by = {Pt, Zt, at}, with P and Z denoting concen- 
trations of phytoplankton and zooplankton, respectively, and a the stochastic growth 
rate of phytoplankton. These interact via: 

atP - cPZ (24) 
ecPZ — miZ — rrigZ'^ . (25) 



dP 

~dt 
dZ 

~dt 
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Here, t is time in days, with prescribed constants c = .25, e = .3, = .1 and = .1. 
The stochastic growth term, at, is drawn daily as at ~ N'{fi,a). Parameters to be 
estimated are © = {/i, a}. Uniform prior distributions are assigned to the parameters, 
/i ~ W(0, 1) and a ~ W(0,.5). Log-normal distributions are placed over the initial 
conditions, InP ~ A/'(ln2, .2) and InZ ~ A/'(ln2, .1). Phytoplankton (P) is observed 
with log-normal noise: 

lnPobs~Ar(lnP,.2). (26) 

Simulated data is used by integrating forward a single trajectory for 100 days, tak- 
ing P daily and adding observation noise. Z is unobserved. 

A joint UKF is first applied, state augmented to include parameters, obtaining a 
preliminary approximation to the posterior over these. This provides the starting dis- 
tribution for Markov chains, as well the covariance (after scaling by 0.18) for a random- 
walk Gaussian proposal distribution. This scaling factor is chosen using pilot runs 
with the PFO method at 64 particles. Starting at the rule-of-thumb 2.4"^ /Ng = 2.88 [9|, it 
is halved (four times) until a mixing rate close to the rule-of-thumb 23% [9] is achieved 
in the first 500 steps of the chain. Initial conditions are importance sampled in the par- 
ticle filter, i.e. the target is factorised according to Mixing might be improved by 
using ([6]); the ulterior motive, already mentioned, is for ready visualisation of Markov 
chain behaviour over only two dimensions. 

Conditional acceptance rates are computed (as described in ^: A I at 1024 points on a 
regular grid across the uniform prior distribution, using 200 likelihood evaluations at 
each point. To produce contours of the surface, a Gaussian process is fit to the points 
by maximum likelihood, with a constant mean function, isotropic squared exponential 
covariance function and Gaussian likelihood [jSTl |30|. Results are presented for the 
particle-matched case in Figure |2| and compute-matched case in Figure |3] 

Multiple chains are run for each method, with the multivariate P*" statistic of Brooks 
and Gelman [2] computed at regular intervals to assess rapidity of convergence (Fig- 
urej4]). The posterior distribution obtained over parameters is marked in Figures |2] 
& 131 using samples drawn across all chains for each method. For a select chain, the 
state posterior is visualised in Figure |5| along with the ground truth trajectory and 
observations for comparison. 



5.2 NPZD model 

The introduction of nutrients, A^, and detritus, D, into the PZ model provides a more 
realistic system with which real observational data can begin to be assimilated: an 
NPZD model. These additional terms are accompanied by various environmental 
forcings and rate processes that produce a more challenging model, with nonlinear 
responses ranging from convergence, to periodicity, to chaos. The full details and 
motivation behind the model are given in Parslow et al. Il25l . A brief description to 
elucidate some of the complexity is given in ^|C] 

A data set is generated by simulating the model with artificial forcing for 100 days, 
from which nutrient (A^) and chlorophyll-a (Chla) observations are produced daily. 

A joint UKF followed by an unscented Rauch-Tung-Striebel smoother [URTSS,l32| 
are applied, taking the joint covariance over initial conditions and parameters (scaled 
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Figure 2: Conditional acceptance rates for each method applied to to the PZ case study, 
particle-matched, across the support of the uniform prior distribution. Darker shading 
denotes higher CAR (see key bottom right). The dots at approximately (0.3, 0.1) in each 
plot mark the ground truth parameters from which the data is simulated, and the bold 
contours nearby the posterior distribution obtained. 
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Figure 3: Conditional acceptance rates for each method applied to to the PZ case study, 
compute-matched. See Figure |2] caption for details. 
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Figure 4: Convergence rates of Markov chains for the PZ case study, (left) particle- 
matched, and (right) compute-matched. Each line shows the evolution of the RP statis- 
tic of Brooks and Gelman [2J for a particular method (see key), as the number of steps 
taken increases. The statistic is computed using 256 chains for each method. A more 
rapid approach to 1 indicates faster convergence. 
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Figure 5: Time-marginal posteriors over state variables (P and Z) and noise term (a) 
for the PZ case study. Results are as obtained by PMMH using the CUPFl method. 
Other methods give comparable results, although not shown. Bold lines of the prior 
and posterior denote the median, and shaded regions the 95% credibility interval, of 
the distribution. 



by .012) for a random-walk Gaussian proposal for Markov chains, and both the joint 
mean and covariance as the starting distribution. Conditional acceptance rates are 
computed (as described in ^: A I at points drawn randomly from the prior distribution. 
Pairwise copulas between parameters and the CAR, computed empirically, are shown 
in Figure |6| with the empirical cdfs of those CAR in Figure |7j Figure |8] shows the 
covariance matrix over parameters obtained by both the PMMH and approximating 
URTSS method. The BP statistic [2J is computed across multiple chains and shown in 
Figure |9] The univariate posterior marginal distributions of parameters and state are 



given in Figures 10 & 11 



6 Discussion 

The collapse of a state-space model to its independent noise terms yields a closed-form 
transition density regardless of the functional forms of the process under study. Once 
this is acquired, the proposal distributions within an auxiliary particle filter (APF) can 
be tuned to improve sampling efficacy. It should be expected that the posterior den- 
sities over such noise terms may exhibit significant skew and kurtosis, for which it is 
difficult to devise good proposal distributions. Nevertheless, for the realistic marine 
biogeochemical models explored in this work, there is demonstrable utility in tuning 
merely the mean and covariance of a multivariate Gaussian proposal, increasing con- 
ditional acceptance rates, reducing variance in the likelihood estimator and improving 
overall convergence rates. 

A common ill of the particle marginal Metropolis-Hastings (PMMH) sampler is 
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Figure 6: Estimated copula functions between parameters (x-axes) and CAR (y-axes) 
for the NPZD case study, using the CUPFl method (other methods give similar re- 
sults). For parameters, the prior univariate cumulative density functions (cdf) are used 
for transformation to uniform marginals; for the CAR, the empirical cdf is used. The 
copula function is approximated using a kernel density estimate of bandwidth .075 
over 4096 points sampled from the prior distribution over parameters, with CARs 
computed from 200 likelihood evaluations at each point. Edge effects are an unfortu- 
nate artifact of the kernel density estimate. Most striking is the significant sensitivity 
of the CAR to the jjiciz parameter, explained in the text. 
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Figure 7: Empirical cdfs of the CAR for each method on the NPZD case study, (left) 
particle-matched, and (right) compute-matched. For each method, the CAR is com- 
puted at the same sample points used to construct Figure [6| then the empirical cdf 
over all of these CARs computed. As higher CARs are preferred, a lower cumulative 
density on the y-axis is preferred for any given rate on the x-axis. 
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Figure 8: Covariance matrices over parameters of the NPZD model, as returned by 
(left) PMMH using CUPFl, and (right) URTSS. The area of each square is propor- 
tional to magnitude, with filled squares denoting positive, and empty squares nega- 
tive, covariance. On the whole, some significant off-diagonal covariance is captured 
by the inexpensive URTSS method, from which a good proposal distribution might be 
constructed. 
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Figure 9: Convergence rates of Markov chains for the NPZD case study, (left) particle- 
matched, and (right) compute-matched. See Figure |4] for explanation. Computed us- 
ing 256 chains for each method. 
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Figure 10: Time-marginal posterior distributions over state variables for the NPZD 
case study. Results are as obtained by PMMH using the CUPFl method. Other meth- 
ods give comparable results, although not shown. Bold lines of the prior and posterior 
denote the median, and shaded regions the 95% credibility interval, of the distribution. 
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for it to enter a state on the basis of an unusually high likelihood estimate, and then 
remain there for a prolonged period. This may manifest a chain that stammers rather 
than mixing fluently The behaviour is more pronounced when the variability of the 
likelihood estimator is high, and further confounded by heteroskedasticity across the 
domain of parameters. This work introduces the conditional acceptance rate (CAR) as a 
useful measure of the impact of such variability on the mixing of a PMMH chain at any 
given point. Applying this at multiple points reveals substantial variability in both the 
PZ (Figures|2]&|3]) and NPZD (Figure|7| cases. For the PZ model, a clear decline in CAR 
away from the posterior mode is apparent (Figures |2] & |3]), although larger values of 
the diffusive parameter a lend improvement. The NPZD case study is more complex, 
many parameters having no apparent influence on CAR (Figure [6]), at least over the 
support of the prior distribution. Again, however, there is some indication that larger 
values of diffusive parameters (specifically ZDF) improve CAR, which also falls away 
dramatically with distance from the posterior mode of a parameter to which the model 
is particularly sensitive, /ic; Jj 




The dependence of CAR on diffusive terms is not surprising given that the pro- 
cess model itself forms part of the APF proposal distribution in all methods explored. 
Broader proposals tend to deliver more effective importance samples, or are at least 
more forgiving of deviation between the modes of the proposal and target distribu- 
tions. Note that the CAR will likely decline again if diffusive parameters become too 
large, but the prior distribution in both the PZ and NPZD cases constrains support 
to regions where a simple correlation is apparent. Because the prior distribution over 
parameters does not factor into any calculations made by the particle filter, the in- 
terpretation of the other apparent dependency of CAR should be distance from the 
maximum likelihood estimate (MLE), rather than distance from the posterior mode, 
or maximum a posterior (MAP) estimate. 

Given the variability of the CAR, there are two potential pitfalls to avoid when 
using a PMMH sampler: 

(a) initialisation in a region where CAR is low, and 

(b) a particularly informative prior distribution that biases the posterior sufficiently 
far from the MLE to be in a region where CAR is low. 

Either case may result in a failure of the Markov chain to converge, or even accept, 
in reasonable time. These should be considered failure modes of the PMMH sam- 
pler in much the same way as strongly correlated variables can cause slow mixing in 
the Gibbs sampler, or multiple modes can cause quasiergodicity in any MCMC proce- 
dure. The CAR is one avenue to diagnose this behaviour. Note that one cannot simply 
reduce the size of the proposal to arbitrarily increase acceptance rates, as with conven- 
tional MCMC with exact likelihoods. The construction of the CAR is as an estimate of 
the acceptance rate of repeated Dirac proposals at a single point of interest: the lower 

^This parameter dictates the mean of the stationary distribution of the zooplankton clearance rate au- 
toregressive. At low clearance rates, phytoplankton (P) will periodically escape zooplankton (Z) graz- 
ing control and begin a rapid bloom, triggering spikes in chlorophyll-a {Chla) that cannot be reconciled 
with observations. At high clearance rates, phytoplankton is relentlessly suppressed by zooplankton 
predation, keeping chlorophyll-a at much lower values than those observed. 
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limit of proposal size. Neither is relying on an arbitrary initialisation, and expecting a 
PMMH chain to converge to the target distribution, likely to succeed: even a smooth 
likelihood surface will appear noisy to this sampler. The aim should be to constrain 
variability to a tolerable level within the support of the posterior distribution, and use 
some approximating method to initialise the Markov chain within sight of this. The 
use of the UKF and URTSS for this purpose appears successful for the case studies 
here. As an aside, the approximate posterior covariance output from these methods 
seems to capture some significant off-diagonal elements, so as to make a more useful 
proposal covariance than a contrived orthogonal structure (Figure |8]). 

When the number of particles (M) is matched across methods, the MUFF and 
CUFF methods outperform the basic PF methods in both CAR (Figure |2| left-hand 
side of Figure [7]) and convergence rate (left-hand side of Figures |4] & |9]). We stress 
again, however, that the process model must mix sufficiently quickly for the proposal 
over noise terms to be of any benefit. This is particularly the case for MUFF, which ne- 
glects information on the current state of each particle. For slow-mixing models with 
long memory, the look-behind is far more important, and the tuning of proposals over 
noise terms may be additional computation for naught. Conversely, the look-behind 
is precisely only of benefit for such slow-mixing models with some memory. Relative 
to the PZ model, the NPZD model fits this description. For fast-mixing models the 
look-behind is of little benefit, and even harmful if the single-point pilot does not rep- 
resent the whole conditional transition density well. This appears to be the case with 
the PZ model, where PFl and MUPFl exhibit smaller CAR than their PFO and MUPFO 
partners, exaggerated as the pilot becomes less representative at higher o. Note that 
CUPFl implements its lookahead using UKFs conditioned on each particle, not using 
single-point pilots, and this should rarely harm results. Indeed, CUPFl requires com- 
parably little additional computation over CUPFO, and is almost certainly worth using 
over it for additional robustness. 

In compute-matched scenarios, the MUFF methods appear to retain some superi- 
ority (Figure |3| right-hand side of Figures |4| [7| & |9]). The CUFF methods do not. We 
should hesitate to dismiss the conditional approach, however, noting that the compu- 
tational burden is from repeated UKFs, which could be replaced with, say, extended 
Kalman filters (EKFs) to deliver similar returns for less computational investment. An 
interesting subplot arises from the CUFF methods in the PZ case study. No other 
method matches the CAR of CUPFO and CUPFl far from the posterior distribution 
(Figure |2]), even after correction for compute time (Figure |3]). This suggests that these 
methods may make a more robust choice for early steps in a Markov chain if a good 
initialisation is unavailable, later displaced by one of the cheaper methods once the 
posterior mode is in sight. Note, finally, that the overall acceptance rates achieved by 
each method are commensurate with their corresponding CARs in the vicinity of the 
posterior mode, rather than CARs across the prior as a whole (most clear in the PZ 
case, compare Table |2] against Figures |2] & |3]). This is, of course, where a successful 
chain will spend the bulk of its time. 

Finally, we return to practical matters on the marine biogeochemical case studies. 
At least one of the operational requirements of such models is their skill in prediction. 
Such a use typically requires slow-mixing process models, not the fast-mixing models 
desired for successful application of the MUFF and CUFF methods presented here. 
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We note, however, that proposals over noise terms provide some benefit for all but 
the slowest mixing of models. It is also worth noting that parameter estimation alone 
can significantly improve long-term predictions over those of the prior model, even 
for fast mixing models. The handling of sparse data is also an issue. With some exten- 
sion, both the MUFF and CUFF methods could be made useful in this context, such 
as for the observational data sets available at the site of OSF. This might involve data 
augmentation II35I , interleaving simulation of missing data values given the current 
state of the chain with a FMMH step conditioned on both real and imputed observa- 
tions. This would mean that the distance to which the UKF (or EKF, or some other 
such method) looks ahead to the next observation (real or imputed) could be limited 
to an interval of time for which the approximation is credible. We have left this to 
future work. 

The collapse of a state-space model to its noise terms seems a generally good ap- 
proach to dealing with the absence of a closed-form transition density, which would 
otherwise prohibit the use of clever proposal distributions in inference. This work 
has demonstrated the utility of doing so to improve the acceptance and convergence 
rates of a FMMH sampler, working the UKF into the proposal mechanism of the in- 
ner particle filter. In doing so it has presented empirical results that elucidate some 
of the behaviours peculiar to the FMMH sampler, particularly highlighting the het- 
eroskedasticity of the likelihood estimator, and the criticality of a good initialisation 
for the successful convergence of the Markov chain. 

A Approximating the conditional acceptance rate (CAR) 

At a particular point z = {xq, 6}, we wish to compute the long-term acceptance rate of 
a particle marginal Metropolis-Hastings (FMMH) chain starting at z, and continuously 
proposing an identity move to z. A finite number, L, of likelihood estimates, /i, . . . , II, 
are made at z using some particle filtering method. From these a transition matrix 
T G M-^^^is formed: 



In the extended space of the FMMH sampler, Tij gives the probability of accepting a 
move to the jth sampled state, of likelihood Ij, from the ith sampled state, of likelihood 
li, under a uniform prior and proposal. From state i, the probability of accepting a 
uniformly drawn proposal is thus 



The 1/L term has the effect of counting a proposal to remain at state i as an acceptance. 

From some arbitrary state, the Markov model defined by T may be run to equilib- 
rium, where the probability of being in state i is simply the normalised term 




ifi=j. 



(27) 



1-Tu + 1/L . 



(28) 



k 



(29) 
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The long-term average acceptance rate, a, is then: 

a =T-/3. 



(30) 



We call this the conditional acceptance rate at z, abbreviated CAR. 

The statistic may be computed more readily when the normalised likelihoods, 
li, . . . Jl, are pre-sorted in ascending order. In this case: 



1 

L 



(31) 



Proceeding from (30 1, let c be the vector of inclusive prefix-sums over 1 (i.e. Cj 
Yl]=i h)' proceed: 



La 



r=l Lj=1 J 
L i L 

i=l j=l i=l 
L L L 

1=1 i=l i=l 

L 

i=l 



(32) 
(33) 
(34) 
(35) 



Thus, each a is readily computed by sorting and normalising 1 to produce I, prefix- 
summing over this to compute c, and summing over this, combined with some simple 
scalar operations, to close. 



B Computational matters 

B.l Asynchronity in the marginal unscented particle filter 

For the MUPF methods, observe that while the PF estimate of ]9(ut|ui;f_i, xq, 0) de- 
pends on prior evaluation of the UKF estimate pj^{ut\ui;t-i, xq, 6) to form a proposal 
distribution, there is no reciprocal dependence of the UKF on the PF. Consequently, 
the UKF may be run asynchronously to the PF, as long as it remains at least one step 
ahead at all times. This is not the case for CUPF, which requires the samples u}:^j^ 
from the PF before forming the proposal distributions for time t. 

This property is particularly advantageous in a heterogeneous hardware context. 
The computational expense of repeated particle filters in PMCMC is offset by the in- 
herent parallelism of the SMC component. The large number of independent parti- 
cle propagations is particularly amenable to implementation on graphics processing 
units [CPUs, HZl |23l |22l, being characterised by SIMD (Single-Instruction Multiple- 
Data) vector operations. The UKF, on the other hand, is dominated by matrix opera- 
tions such as the Cholesky factorisation, for which GPU implementations [20J tend to 
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be favourable over central processing unit (CPU) implementations only for very many 
dimensions [36 J. For a moderate number of dimensions, on the order of a few hun- 
dred, CPU execution of the UKF is preferred. If UKF execution time on the CPU is less 
than PF execution time on the GPU, MUPF may be employed at negligible overhead 
to the ordinary bootstrap particle filter, assuming that the CPU would otherwise idle 
during GPU execution. 



B.2 Weight computations for Gaussian noise and proposal 

Consider the weight computation for the auxiliary particle filter in ( [T0| >, and in par- 
ticular the factor p{ut)/qt{\it), where p{ut) = i.i.d. A/'(0, 1) and qt{ut) = A/'(^i, S). For 
models with Gaussian noise, this case includes the MUPF and CUPF methods intro- 
duced in ^|4l 



To sample u^^ from g2(-) in (11 1, one draws ^™ ~ i.i.d. A/'(0, 1) and sets = L^J" + 



fi, where L denotes the lower-triangular Cholesky factor of E. Let 

\nvr = \np{ut)-\nqt{ut) (36) 

= ^in|si-^(ur)^ur + ^(ur-/x)^s-i(ur-M) 07) 

= \n\L\ + l[iCfC-i<f<]- (38) 



Substituting i;™ for p{ut) / (uj) in ( 10 1 gives the appropriate weighting for particle m. 

In the case of the MUPF, where gt(ut) = p_^r(uj|xo, 0, yi;*) = J^{n = = Si), 

continued propagation of the UKF requires computation of the Cholesky factor L = 
Lt anyway, so that the additional operations required for the weight adjustment are 
comparably negligible: two dot products and a single vector multiply-and-add for 
each particle, and a single sum over the diagonal of Lt to compute In |Lt|. These may, 
of course, be combined into matrix operations to perform the computations for all 
particles together. 

For the CUPF, denoting k = a!^ for brevity, gt(u]") = pAr(u™|u5'.i_i, xq, 0, yt) = 

N'{^i = Uf, S = S^), and each UKF need only deliver fi^ and after conditioning 
on yt. Producing is particularly cheap by exploiting the fact that the uncorrected 

prediction of is always I, giving the measurement update: 



^ -C^,(S^)-^(C^,)^ (39) 

^k /J k\ — l /J k\—T f f~ik 



^-cl,{iJ;)-\l^)-\ciy, (40) 



" k 

where = cov(Yt, Yt|u^.j_]^, xq, 0), and C^^^ = cov(Ui, Yt|u^.^_]^, xq, 0). Because the 

" k 

Cholesky factor of I is of course I, that of the updated Sj may be computed using Ny 
rank-1 Cholesky downdates to I over the columns of (L^)~^(C^y)^. For small Ny, this 
will be substantially quicker than the full computation and Cholesky factorisation of 

" k " k 

the updated . It may be possible to further exploit Sj = I in this operation, such as 
using sparse formulations BH, but we are yet to explore the possibilities thoroughly. 
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C NPZD model specification 



The NPZD model represents the interaction of nutrients (N), phytoplankton (P), zoo- 
plankton (Z) and detritus (D), in the common currency of nitrogen, within the surface 
mixed layer of a body of water. The surface waters are modelled as a single box, 
subject to exogenous environmental forcings such as available light, temperature and 
changes in mixed layer depth. A full derivation of the model, including sources for the 
prior distribution over parameters, can be found in Parslow et al. [25J; the essentials 
for the present work are summarised here. We have adopted a slightly simpler version 
of the model than that presented in Parslow et al. II25L whereby exogenous forcings 
are considered time-invariant. 

For reference as components of the model are introduced, the model state is X = 
{P, Z, D, N, g'"'''', A"^^^ Piv, aN, Iz, Clz, Ez, r^, tuq}, and the parameters = {Kw, 

O'Ch-, Sd, fn, /igmax, fix'^^'', ^J'Rn, f^arf: fJ'Iz^ f^Clzy f'Ez: f^rny f^mg, PDF, ZDF}. 

C.l Noise model 

The model features nine noise terms, for i = 1, . . . , 9, each coupled to a univariate 
autoregressive process Pj. Four of these are phytoplankton-related, given by 

Pi(t + At) = P,(t) ■ (1 - At/rp) + ifii + PDF ■ aS) ■ At/rp , (41) 

where At is a discrete time step (one day), fii a parameter to be estimated, PDF a com- 
mon diversity factor to be estimated, (Xj a prescribed scaling factor, and Tp a common 
characteristic time scale, also prescribed. The remaining five autoregressive processes 
are zooplankton-related, modelled using the same form, with ZDF and tz replacing 
PDF and rp, respectively. 

Each process represents a property of the phytoplankton (zooplankton) commu- 
nity, the species composition of which will change with time. Rather than model indi- 
vidual species, the phytoplankton (zooplankton) community is modelled collectively, 
with diversity factors PDF and ZDF scaling stochastic drivers used to model the 
changing influence of community composition. 

The four phytoplankton processes are {g'^^^, A'"^", Rn, aj^}, and the five zooplank- 
ton processes {Iz,Clz, Ez^rp), mq}. Each is accompanied by its matching noise term 
amongst U = {^^max, ^;,max, , , ^j^, ^ciz, ^Ez, ^td, ^mg} , and mean parameter 

amongst {/Xgmax, /i^max, fj^p^, /Xq^, /i/^ , flQlz^ f^Ez^ /^r^ , /^mg}- 
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C.2 State model 



The equations governing interactions between the remaining state variables {P, Z, D, N} 



Ez- gr- Z -m- Z (43) 



are: 

dP K 

H = 0-P-,r.Z^ — .iBCP-P) (42) 

dZ 
~dt 

f = (l-E,).fo-<,r.Z + r„.Z-r.D-So-j^ + j^-(BCD-mi) 

dN K 

— = -g.p + {l-Ez)-{l-fD)-9r-Z + r-D + j^-iBCN-N). (45) 

Here, g is the phytoplankton specific growth rate (per day, or d"^), gr is the zooplank- 
ton specific grazing rate (mg P grazed per mg Z d~^), m is the zooplankton specific 
mortality rate (d^), and r is the specific breakdown rate of detritus (d^). The fraction 
Ez (a parameter to be estimated) of zooplankton ingestion is converted to zooplankton 
growth and, of the remainder, a fraction /d allocated to detritus and the rest released 
as dissolved inorganic nutrient, A^. The grazing rate gr, mortality rate m and growth 
rate g are not only functions of the state, but also prescribed exogenous forcings and 
physiological constants, described below. 

C.3 Rate processes 

A multiplicative temperature correction Tc is applied to all rate processes, for which a 
Qio formulation for dependence on temperature, T, is used: 

Tc = gSJ-^-^^/'°, (46) 

where T^f is a reference temperature, and Qio a prescribed constant. 

The zooplankton grazing rate (gr) is dependent on the phytoplankton concentra- 
tion (zooplankton functional response): 

gr = ^'-'^-^^ (47) 

where v is a given power; the relative availability of phytoplankton. A, is 

A=^; (48) 

Iz is the maximum zooplankton ingestion rate (mg P per mg Z per day); Clz is the 
maximum clearance rate (volume in m"^ swept clear per mg Z per day). For v = 1, (47 1 



takes the form of a Type-2 functional response (standard rectangular hyperbola) [13.1 , 
and for v > 1 a Type-3 sigmoid functional response. 

A quadratic formulation for zooplankton mortality is adopted after Steele II33II and 
Steele and Henderson ||34|: 

m = Tc ■ niQ ■ Z , (49) 
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where the quadratic mortality rate uiq has units of d ^(mgZm ^) ^. The detrital rem- 
ineralization rate is dependent only on temperature: 

r = Tc ■ r^i , (50) 

where td prescribes the remineralisation rate at a reference temperature. 

The phytoplankton specific growth rate, g, depends on temperature, T, available 
light or irradiance, E, and dissolved inorganic nutrient, N. It is expressed in terms 
of a maximum specific growth rate at the reference temperature, g™^^ (d^), a light- 
limitation factor. He, and a nutrient-limitation factor, h^'- 

g = Tc-g'^'^-hE-hN/{hE + hN). (51) 
The light-limitation factor is given by 

/lij = 1 - exp(-a ■ A""^" • ^/^'""') , (52) 

where a is the initial slope of the photosynthesis versus irradiance curve (mg C mg 
Chla^^ mol photon^^ m^), and A"^^" is the maximum Chla : C ratio (mg Chla mg C^^). 
Here, a is calculated as the product of the chlorophyll-specific absorption coefficient 
for phytoplankton, ach (m^ mg Chla^^), and the maximum quantum yield for photo- 
synthesis, Q (mg C mol photons^^). E is the mean photosynthetic available radiation 
(PAR) in the mixed layer and is given by 

E = Eo.{l - exp{-Kz))/Kz, (53) 

where Eq is the mean daily photosynthetically available radiation (PAR) just below the 
air-sea interface and Kz is given by 

Kz = {Kw + ach ■ Chla) ■ MLD. (54) 

In ( [54| , Kw is attenuation due to the seawater and ach- 
The nutrient-limitation factor is given by 

= (^?-- ■ TcM + N ' ^^^^ 

where aj^ is the maximum specific affinity for nitrogen uptake (d~^ mg A^~^ m^). 
The phytoplankton N : C ratio, x, predicted by the model is given by 

riE + riN 

where x"^'^ and x'"^" are the prescribed minimum and maximum : C ratios (mg 
mg C-^). 



C.4 Boundary conditions 

The simple one-box mixed layer model adopted here needs to allow for the effects of 
physical exchanges between the mixed layer and the underlying water mass. With the 
exception of BCN , all boundary conditions {BCP, BCD, BCZ) are set to zero for the 
experiments in this work. The variable k sets the strength of the mixing; in this study, 
we assume that the lower two metres of the mixed layer are replenished daily with 
water from the sub-mixed layer. MLD and BCN are in this case time invariant and 
set to 40 m and 200 mg N m'^ respectively. 
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C.5 Observation model 



The model predicts the phytoplankton Chla: C ratio A, and this can be combined with 
the N : C ratio x to convert phytoplankton biomass P (mg N m~^) to a predicted Chla 
concentration: 

Chla = P ■ (A'^^Vx'""") ■ hN ■ Tc/{Rn ■ hs + /ijv) • (57) 

Both N and Chla are observed, each with log-normal noise of 40%, i.e. IniVobs ^ 
A/'(ln A^, .4), and In Chlaobs ~ A/'(ln Chla, A). Observations can thus be written as Y = 

{iVobs, Chlaobs}- 
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